// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>

#include <cmdioutil.h>
#include "xyz.h"


void get_index_list(const util::list_expression<int> &expr, int noa,std::list<int>& index_list)
{
        for(util::list_expression<int>::const_iterator it = expr.begin(); it != expr.end(); it++)
        if(*it == 0) throw std::runtime_error("Index cannot be zero in the index list.");
        else if((*it > noa) || (*it < -noa)) throw std::runtime_error("Index is out of range in the index list. There are less atoms!");
        else
        {
            int index = *it;
            int i = index < 0 ? index + noa : index - 1;
            index_list.push_back(i);
        }
}

void get_index_list(int ia,int na,std::list<int>& index_list)
{
        index_list.clear();
        for(int i = ia; i < ia+na; i++) 
        {
            index_list.push_back(i);
        }
}

struct SH
{
    Eigen::MatrixXcd S,H;
    double ef;
    void resize(int r,int c) { S.resize(r,c); H.resize(r,c); }
};

struct blockSH : SH
{
    std::vector<std::pair<int,int> > orbital_range;
    int n;
    void get_block(int ii,int nii,int jj,int njj,SH & sh) const
    {
        ii = ii < 0 ? ii + n : ii - 1;
        jj = jj < 0 ? jj + n : jj - 1;
        
    
        int nI,nJ;
        nI = 0; for(int i = ii; i < ii+nii; i++) nI += orbital_range[i].second;
        nJ = 0; for(int j = jj; j < jj+njj; j++) nJ += orbital_range[j].second;

        sh.resize(nI,nJ);
        int li = 0;
        for(int I = ii; I < ii+nii; I++)
        {
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(int J = jj; J < jj+njj; J++)
            {
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                sh.S.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                sh.H.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
        sh.ef = ef;
    }
        
    void get_part(int ii,int nii,SH & h) const
    {
        ii = ii < 0 ? ii + n : ii - 1;


        int nI;
        nI = 0;
        for(int i = ii; i < ii+nii; i++) nI += orbital_range[i].second;

        sh.resize(nI,nI);
        int li = 0;
        for(int I = ii; I < ii+nii; I++)
        {
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(int J = ii; J < ii+nii; J++)
            {
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                sh.S.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                sh.H.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
        sh.ef = ef;
    }

    void get_part(int ii1,int nii1,int ii2,int nii2,SH & sh11,SH & sh12,SH & sh21,SH & sh22) const
    {
        get_part(ii1,nii1,sh11);
        get_part(ii2,nii2,sh22);
        get_block(ii1,nii1,ii2,nii2,sh12);
        get_block(ii2,nii2,ii1,nii1,sh21);
    }
        
    void get_atom(int I,int J,SH & sh) const
    {
        I = I < 0 ? I + n : I - 1;
        J = J < 0 ? J + n : J - 1;

        int i,j,ni,nj;


        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        j = orbital_range[J].first;
        nj = orbital_range[J].second;
        sh.S = S.block(i,j,ni,nj);
        sh.H = H.block(i,j,ni,nj);
        
        sh.ef = ef;
    }

    void get_block(const std::list<int> & index_list_I,const std::list<int> & index_list_J,SH & sh) const
    {
        int nI,nJ;
        nI = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) nI += orbital_range[*it].second;
        nJ = 0;
        for(std::list<int>::const_iterator it = index_list_J.begin(); it != index_list_J.end(); it++) nJ += orbital_range[*it].second;
    
        sh.resize(nI,nJ);
        int li = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
        {
            int I = *it;
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(std::list<int>::const_iterator jt = index_list_J.begin(); jt != index_list_J.end(); jt++) 
            {
                int J = *jt;
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                sh.S.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                sh.H.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    }
    
    
    
/*    
    void set_orbital_range(int natom,int norb)
    {
        orbital_range.clear();
        for(int i=0; i<natom;i++)
        orbital_range.push_back(std::pair<int,int>(i*norb,norb));
        S.resize(natom*norb,natom*norb);    
        H.resize(natom*norb,natom*norb);    
        ef=0.0;
    }
*/
    
    void get_block(int I,int J,Eigen::MatrixXcd & s,Eigen::MatrixXcd & h) const
    {
        int i,j,ni,nj;
        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        j = orbital_range[J].first;
        nj = orbital_range[J].second;
        s = S.block(i,j,ni,nj);
        h = H.block(i,j,ni,nj);
    }

    void set_block(int I,int J,const Eigen::MatrixXcd & s,const Eigen::MatrixXcd & h)
    {
        int i,j,ni,nj;
        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        j = orbital_range[J].first;
        nj = orbital_range[J].second;
        S.block(i,j,ni,nj) = s;
        H.block(i,j,ni,nj) = h;
    }

    void get_block(const std::list<int> & index_list_I,const std::list<int> & index_list_J,Eigen::MatrixXcd & s,Eigen::MatrixXcd & h) const
    {
        int nI,nJ;
        nI = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) nI += orbital_range[*it].second;
        nJ = 0;
        for(std::list<int>::const_iterator it = index_list_J.begin(); it != index_list_J.end(); it++) nJ += orbital_range[*it].second;
    
        s.resize(nI,nJ);
        h.resize(nI,nJ);
        int li = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
        {
            int I = *it;
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(std::list<int>::const_iterator jt = index_list_J.begin(); jt != index_list_J.end(); jt++) 
            {
                int J = *jt;
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                s.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                h.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    }




    void get_partition(const std::list<int> & index_list_I,blockSH & sh) const
    {
        int nI;
        nI = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) nI += orbital_range[*it].second;
    
        sh.orbital_range.clear();
        sh.S.resize(nI,nI);
        sh.H.resize(nI,nI);
        int li = 0;
        int inc = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
        {
            int I = *it;
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            sh.orbital_range.push_back(std::pair<int,int>(inc,ni)); inc+=ni;

            int lj = 0;
            for(std::list<int>::const_iterator jt = index_list_I.begin(); jt != index_list_I.end(); jt++) 
            {
                int J = *jt;
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                sh.S.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                sh.H.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    }
    

    void get_partition(int ii,int nii,blockSH & sh) const
    {
        int nI;
        nI = 0;
        for(int i = ii; i < ii+nii; i++) nI += orbital_range[i].second;

        sh.orbital_range.clear();
        sh.S.resize(nI,nI);
        sh.H.resize(nI,nI);
        int li = 0;
        int inc = 0;
        for(int I = ii; I < ii+nii; I++)
        {
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            sh.orbital_range.push_back(std::pair<int,int>(inc,ni)); inc+=ni;

            int lj = 0;
            for(int J = ii; J < ii+nii; J++)
            {
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                sh.S.block(li,lj,ni,nj) = S.block(i,j,ni,nj);
                sh.H.block(li,lj,ni,nj) = H.block(i,j,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    
    }

    
    void add_partition(blockSH & sh) 
    {
        int nI;
        nI = 0;
        for(std::vector<std::pair<int,int> >::iterator it = orbital_range.begin(); it != orbital_range.end(); it++) nI += it->second;
        
        for(std::vector<std::pair<int,int> >::iterator it = sh.orbital_range.begin(); it != sh.orbital_range.end(); it++)
        {
            int i,ni;
            i = it->first;
            ni = it->second;

            orbital_range.push_back(std::pair<int,int>(i+nI,ni));
        }
        Eigen::MatrixXcd ss,hh;
        ss.resize(S.rows()+sh.S.rows(),S.cols()+sh.S.cols());
        hh.resize(H.rows()+sh.H.rows(),H.cols()+sh.H.cols());
        ss.setZero();
        hh.setZero();
        
        ss.block(0,0,S.rows(),S.cols()) = S;
        hh.block(0,0,H.rows(),H.cols()) = H;

        ss.block(S.rows(),S.cols(),sh.S.rows(),sh.S.cols()) = sh.S;
        hh.block(H.rows(),H.cols(),sh.H.rows(),sh.H.cols()) = sh.H;
        
        S = ss;
        H = hh;        
    }


    void set_block(const std::list<int> & index_list_I,const std::list<int> & index_list_J,const Eigen::MatrixXcd & s,const Eigen::MatrixXcd & h)
    {
        int li = 0;
        for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
        {
            int I = *it;
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(std::list<int>::const_iterator jt = index_list_J.begin(); jt != index_list_J.end(); jt++) 
            {
                int J = *jt;
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                S.block(i,j,ni,nj) = s.block(li,lj,ni,nj);
                H.block(i,j,ni,nj) = h.block(li,lj,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    }



    void set_block(int ii,int nii,int jj,int njj,const Eigen::MatrixXcd & s,const Eigen::MatrixXcd & h)
    {
        int li = 0;
        for(int I = ii; I < ii+nii; I++)
        {
            int i,ni;
            i = orbital_range[I].first;
            ni = orbital_range[I].second;
            int lj = 0;
            for(int J = jj; J < jj+njj; J++)
            {
                int j,nj;
                j = orbital_range[J].first;
                nj = orbital_range[J].second;
                S.block(i,j,ni,nj) = s.block(li,lj,ni,nj);
                H.block(i,j,ni,nj) = h.block(li,lj,ni,nj);
                lj += nj;
            }
            li += ni;
        }
    }



    void print_pattern_H(std::ostream & os,double tol = 1e-6) const
    {
        for(std::vector<std::pair<int,int> >::const_iterator it = orbital_range.begin(); it != orbital_range.end(); it++) 
        {
            for(std::vector<std::pair<int,int> >::const_iterator jt = orbital_range.begin(); jt != orbital_range.end(); jt++)
            {
                int i,j,ni,nj;
                i = it->first;
                ni = it->second;
                j = jt->first;
                nj = jt->second;
                if(H.block(i,j,ni,nj).norm() > tol)  os << 'X'; else os << '.';
            }
            os << std::endl;
        }
    }
};

void load(std::string fn,blockSH & sh)
{
    std::complex<double> ii(0,1.0);
    
    std::vector<std::string> lines;
    io::load(fn,"#",lines);

    sh.ef = 0.0;
    std::vector<std::string> efermi_block;
    io::parse_block("<fermi-energy>","</fermi-energy>",lines,efermi_block);
    for(std::vector<std::string>::iterator it = efermi_block.begin(); it != efermi_block.end(); it++)
    {
        std::istringstream iss(*it);
        if (!(iss >> sh.ef)) throw std::runtime_error("Not able to parse fermi energy block in the input file!");
        break;
    }

    std::vector<std::string> orbital_block;
    if(!io::parse_block("<index-orbital>","</index-orbital>",lines,orbital_block)) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    if(!orbital_block.size()) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    
    sh.orbital_range.clear();
    int index,orbital;
    int inc = 0;
    int orbital_index = 0;
    for(std::vector<std::string>::iterator it = orbital_block.begin(); it != orbital_block.end(); it++)
    if(util::match("% %",*it,index,orbital))
    {
            std::pair<int,int> p(orbital_index,orbital);
            orbital_index += orbital;
            sh.orbital_range.push_back(p);
            if(++inc != index) throw std::runtime_error("Index has to be in order from 1 in the index-orbital block!");
    }
    int n = orbital_index;
    sh.H.resize(n,n); sh.H.setZero(); sh.S = sh.H; sh.S.setIdentity();
    double re,im;
    int i,j;

    std::vector<std::string> overlap_block;
    if(!io::parse_block("<overlap>","</overlap>",lines,overlap_block)) throw std::runtime_error("There is no overlap matrix block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = overlap_block.begin(); it != overlap_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.S(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!");  }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.S(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!"); }
    
    std::vector<std::string> hamiltonian_block;
    if(!io::parse_block("<hamiltonian>","</hamiltonian>",lines,hamiltonian_block)) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = hamiltonian_block.begin(); it != hamiltonian_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.H(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) sh.H(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    if(!hamiltonian_block.size()) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");

}

void save(std::string fn,blockSH & sh)
{
    std::complex<double> ii(0,1.0);
    std::ofstream os(fn.c_str());
    
    os << "<fermi-energy>" << std::endl;
    os << sh.ef << std::endl;
    os << "</fermi-energy>" << std::endl;
    os << std::endl;
    os << "<index-orbital>" << std::endl;
    for(int i=0; i<sh.orbital_range.size(); i++)    
    os << i+1 << " "  << sh.orbital_range[i].second << std::endl;
    os << "</index-orbital>" << std::endl;
    os << std::endl;
    os << "<overlap>" << std::endl;
    for(int i=0;i<sh.S.rows();i++)
    for(int j=0;j<sh.S.cols();j++)
    if(std::abs(sh.S(i,j))> 1e-7 )
    {
        os << "   " << i+1 << "    " << j+1 << "     "<< std::scientific << std::real(sh.S(i,j)) << "    " << std::scientific << std::imag(sh.S(i,j)) << std::endl;
    }
    os << "</overlap>" << std::endl;
    os << std::endl;
    os << "<hamiltonian>" << std::endl;
    for(int i=0;i<sh.H.rows();i++)
    for(int j=0;j<sh.H.cols();j++)
    if(std::abs(sh.H(i,j))> 1e-7 )
    {
        os << "   " << i+1 << "    " << j+1 << "     "<< std::scientific << std::real(sh.H(i,j)) << "    " << std::scientific << std::imag(sh.H(i,j)) << std::endl;
    }
    os << "</hamiltonian>" << std::endl;

    os.close();
    
}


void load_gollum(std::string fn,blockSH & sh)
{
    std::ifstream is(fn.c_str());
    std::complex<double> ii(0,1.0);
    std::string dump;
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump >> dump)) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning!");
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump >> sh.ef )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning!");
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax errour near the beginning!");
    int rows,cols; 
    if(!(is >> dump >> dump >> rows )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning");
    if(!(is >> dump >> dump >> cols )) throw std::runtime_error("Wrong gollum input format: Syntax error near the beginning");
    int index;
    int pindex = 1;
    int pi = 0;
    int n = rows;    
    sh.orbital_range.clear();
    for(int j = 0; j < n; j++)
    {
        if(!(is >> index >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax error in the orbital block!");
        if(index != pindex)
        {
            std::pair<int,int> p(pi,j-pi);
            sh.orbital_range.push_back(p);
            pindex = index;
            pi = j;
        }
    }
    std::pair<int,int> p(pi,n-pi);
    sh.orbital_range.push_back(p);

    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    //int rows,cols;
    if(!(is >> dump >> dump >> rows  )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    if(!(is >> dump >> dump >> cols  )) throw std::runtime_error("Wrong gollum input format: Syntax error near k-block!");
    for(int j = 0; j < rows; j++)
        if(!(is >> dump >> dump >> dump  )) throw std::runtime_error("Wrong gollum input format: Syntax error in k-block!");
    if(rows > 1) std::cerr << "Warning: Gollum input delivers multiple k point matrix. This load does not handle that, instead try separate xml inputs for each k-point." << std::endl;
    if(!(is >> dump >> dump >> dump >> dump >> dump >> dump  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    //int rows,cols;
    if(!(is >> dump >> dump >> rows  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    if(!(is >> dump >> dump >> cols  )) throw std::runtime_error("Wrong gollum input format: Syntax error near matrix block!");
    sh.H.resize(n,n); sh.H.setZero(); sh.S = sh.H;
    for(int k = 0; k < rows; k++)
    {
        int i,j;
        double res,ims,reh,imh;
        if(!(is >> dump >> i >> j >> res >> ims >> reh >> imh  )) throw std::runtime_error("Wrong gollum input format: Syntax error in matrix block!");
        sh.H(i-1,j-1)=reh+ii*imh;
        sh.S(i-1,j-1)=res+ii*ims;

    }
    is.close();
    
}



void get_cmd(int argc, char* argv[],blockSH & sh,std::string & out)
{
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::posarg<std::string> hamla_fn(arg,"input a-lead hamiltonian file name", "");
        cmd::posarg<int> aoa(arg,"exmol lead atom index overlap", 0);
        cmd::posarg<std::string> hamm_fn(arg,"input exmol hamiltonian file name", "");
        cmd::posarg<int> aob(arg,"exmol lead atom index overlap", 0);
        cmd::posarg<std::string> hamlb_fn(arg,"input b-lead hamiltonian file name", "");
        cmd::posarg<std::string> out_fn(arg,"out hamiltonian file name", "");
        
        if(*help)
        {
                arg.print_help(std::cout,"shmat");
                std::cout << "examples:" << std::endl;
                std::cout << "version 0.9b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.9b" << std::endl;
                std::cout << "Copyright (c) 2010 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
        }                                                                                                                


        if(!~hamla_fn) throw std::runtime_error("Give a-lead hamiltonian file!");
        if(!~hamlb_fn) throw std::runtime_error("Give b-lead hamiltonian file!");
        if(!~hamm_fn) throw std::runtime_error("Give extended scatterer hamiltonian file!");
        
        if(!~aoa) throw std::runtime_error("Give number of overlapping atoms with lead a!");
        if(!~aob) throw std::runtime_error("Give number of overlapping atoms with lead b!");

//        std::cout << "ok-1" << std::endl;

        blockSH shla,shlb,shm;
        int ao = *aoa;
        int bo = *aob;
        load(*hamla_fn,shla);
        load(*hamm_fn,shm);
        load(*hamlb_fn,shlb);
    
//        for (int i =0; i< shla.orbital_range.size(); i++)
//        std::cout << shla.orbital_range[i].first << " " << shla.orbital_range[i].second << std::endl;

        
//        blockSH sh;
        int la = shla.orbital_range.size();
        int lm = shm.orbital_range.size();
        int lb = shlb.orbital_range.size();

//        std::cout << "ok0" << std::endl;

        double dea,deb;
        Eigen::MatrixXcd ds,dh,dsm,dhm;
//        std::cout << "la:" << la << "    ao:" << ao << std::endl;    
        
        shla.get_block(la-ao,ao,la-ao,ao,ds,dh);
//        std::cout << "ok0.1" << std::endl;

        shm.get_block(0,ao,0,ao,dsm,dhm);
//        std::cout << "ok0.2" << std::endl;

        dea = std::real((dh-dhm).trace()/ds.trace());
        
//        std::cout << "eshifta : " << dea << std::endl;
        
//        std::cout << "ok0.5" << std::endl;

        shm.get_block(lm-bo,bo,lm-bo,bo,dsm,dhm);
        shlb.get_block(0,bo,0,bo,ds,dh);
        deb = std::real((dh-dhm).trace()/ds.trace());

//        std::cout << "eshiftb : " << deb << std::endl;

//        std::cout << "ok1" << std::endl;

        shla.H = shla.H - dea*shla.S;
        shlb.H = shlb.H - deb*shlb.S;

//        shm.print_pattern_H(std::cout);

        
        blockSH shmm;
        shm.get_block(ao,lm-ao-bo,ao,lm-ao-bo,shmm);
        sh = shla;
//        sh.print_pattern_H(std::cout);
        
//        std::cout << " ---------------- " << std::endl;
//        shmm.print_pattern_H(std::cout);

//        std::cout << "ok2" << std::endl;

        sh.add_block(shmm);

//        sh.print_pattern_H(std::cout);

//        std::cout << "ok3" << std::endl;


        sh.add_block(shlb);

//        std::cout << "ok4" << std::endl;

//        sh.print_pattern_H(std::cout);


        Eigen::MatrixXcd s,h;

        shm.get_block(0,ao,ao,lm-ao-bo,s,h);
//        std::cout << "ok4.5" << std::endl;
//        std::cout << " :lmm: " << lm-ao-bo << std::endl;

        sh.set_block(la-ao,ao,la,lm-ao-bo,s,h);
        
        
//        std::cout << "ok5" << std::endl;

        
        shm.get_block(ao,lm-ao-bo,0,ao,s,h);
        sh.set_block(la,lm-ao-bo,la-ao,ao,s,h);
//        std::cout << "ok6" << std::endl;

        shm.get_block(ao,lm-ao-bo,lm-bo,bo,s,h);
        sh.set_block(la,lm-ao-bo,la+lm-ao-bo,bo,s,h);
//        std::cout << "ok7" << std::endl;

        shm.get_block(lm-bo,bo,ao,lm-ao-bo,s,h);
        sh.set_block(la+lm-ao-bo,bo,la,lm-ao-bo,s,h);

        sh.ef = shm.ef;
        
//        sh.print_pattern_H(std::cout);

        out = *out_fn;        
}


using namespace std;


int main(int argc, char* argv[])
{
try
{
    blockSH sh;
    std::string outfn;
    get_cmd(argc,argv,sh,outfn);
    
    save(outfn,sh);
    
} catch(const char* message)
{
    std::cerr << message << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
catch(const std::exception & e)
{
    std::cerr << e.what() << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
    return EXIT_SUCCESS;



}


/*
int main(int argc, char* argv[])
{
try
{
         rH,rS,H,S;
        
        std::vector<Eigen::Vector3d> rc,c; 
        std::vector<std::pair<int,int> > ror;
        double ref;
        
        get_cmd(argc,argv,rc,ror,rh,rs,ref,c);
        
        Eigen::MatirxXi rh,rk,rl,h,k,l;
        Eigen::Vector3d a1,a2,a3;
        structure_hkl(rc,a1,a2,a3,rh,rk,rl);
        structure_hkl(c,a1,a2,a3,h,k,l);
        
        int rn,n;
        rn = rc.size();
        n = c.size();
    
        //> fast solution
        Eigen::MatrixXcd Hblock,Sblock;
        get_atomic_block_IJ(rH,ror,0,0,Hblock);
        H.resize(n*Hblock.rows(),n*Hblock.cols());
        S.resize(n*Hblock.rows(),n*Hblock.cols());
        H.setZero();
        S.setZero();
        
        //< end
            
        for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        {
            int h_ = h(i,j);
            int k_ = k(i,j);
            int l_ = l(i,j);
            
            bool found = false;
            for(int ri=0; ri<rn; ri++)
            for(int rj=0; rj<rn; rj++)
            if( (h_ == rh(ir,rj))  && (k_ == rk(ir,rj)) && (l_ == rl(ir,rj)) )
            {
                get_atomic_block_IJ(rH,ror,ri,rj,Hblock);
                get_atomic_block_IJ(rS,ror,ri,rj,Sblock);
                found = true;
                break;
            }
            
            if(!found)
            {
                Hblock.setZero();
                Sblock.setZero();
            }
            
            H.block(i*Hblock.rows(),j*Hblock.cols(),Hblock.rows(),Hblock.cols()) = Hblock;
            S.block(i*Hblock.rows(),j*Hblock.cols(),Hblock.rows(),Hblock.cols()) = Sblock;
        
        }

        
} catch(const char* message)
{
    std::cerr << message << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
catch(const std::exception & e)
{
    std::cerr << e.what() << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
    return EXIT_SUCCESS;



}

void coord2nnn(const std::vector<Eigen::Vector3d> & coord,int Iref,int Ia,int Ib,int Ic,std::vector<Eigen::Vector3i> & nnn)
{
    Eigen::Vector3d rref = coord[Iref];
    Eigen::Matrix3d iabc,abc;
    abc.col(0) = coord[Ia]-rref;
    abc.col(1) = coord[Ib]-rref;
    abc.col(2) = coord[Ic]-rref;
    iabc = abc.inverse();

    nnn.clear();
    for(std::vector<Eigen::Vector3d>::const_iterator it = coord.begin(); it != coord.end(); it++)
    {
        Eigen::Vector3i n;
        Eigen::Vector3d r = *it - rref;
        r2n(iabc,r,n);
        nnn.push_back(n);
    }
    
}






void test(const std::vector<Eigen::Vector3d> & coord_model,
          int Iref_model,int Ia_model,int Ib_model,int Ic_model,
          const Eigen::Vector3d & A1,const Eigen::Vector3d & A2,
          const std::list<int> & atom_index_model,
          const std::vector<std::pair<int,int> > & orbital_range_model,
          Eigen::MatrixXcd & H_model,Eigen::MatrixXcd & S_model,
          
          const std::vector<Eigen::Vector3d> & coord_lattice,
          int Iref_lattice,int Ia_lattice,int Ib_lattice,int Ic_lattice,
          const std::vector<std::pair<int,int> > & orbital_range_lattice,
          const Eigen::MatrixXcd & H_lattice,const Eigen::MatrixXcd & S_lattice,
          double rcutoff)
{
    Eigen::Matrix3d iabc_lattice,abc_lattice;
    abc_lattice.col(0) = coord_lattice[Ia_lattice]-coord_lattice[Iref_lattice];
    abc_lattice.col(1) = coord_lattice[Ib_lattice]-coord_lattice[Iref_lattice];
    abc_lattice.col(2) = coord_lattice[Ic_lattice]-coord_lattice[Iref_lattice];
    iabc_lattice = abc_lattice.inverse();
    
    std::vector<Eigen::Vector3i> nnn_lattice;
    for(std::vector<Eigen::Vector3d>::const_iterator it = coord_lattice.begin(); it != coord_lattice.end(); it++)
    {
        Eigen::Vector3d r = *it - coord_lattice[Iref_lattice];
        if(r.norm() < rcutoff)
        {
            Eigen::Vector3i n;
            r2n(iabc_lattice,r,n);
            nnn_lattice.push_back(n);
        }
    }



    Eigen::Matrix3d iabc_model,abc_model;
    abc_model.col(0) = coord_model[Ia_model]-coord_model[Iref_model];;
    abc_model.col(1) = coord_model[Ib_model]-coord_model[Iref_model];;
    abc_model.col(2) = coord_model[Ic_model]-coord_model[Iref_model];;
    iabc_model = abc_model.inverse();
    for(std::list<int>::const_iterator it = atom_index_model.begin(); it != atom_index_model.end(); it++)
    {
        for(double m1 = -1.0; m1 <= 1.0; m1++) 
        for(double m2 = -1.0; m2 <= 1.0; m2++) 
        {
            Eigen::Vector3d s = m1*A1 + m2*A2;
                for(std::list<int>::const_iterator jt = atom_index_model.begin(); jt != atom_index_model.end(); jt++)
            {
                Eigen::Vector3d r = coord_model[*jt]+s - coord_model[*it];
                if(r.norm() < rcutoff)
                {
                    Eigen::Vector3i n;
                    r2n(iabc_model,r,n);
                    int J;
                    if(findnnn(nnn_lattice,n,J))
                    {
//                        H_model[*it,*jt] = H_lattice[Iref_lattice,J];
//                        S_model[*it,*jt] = S_lattice[Iref_lattice,J];
                    } else
                    {
//                        H_model[*it,*jt] = 0;
//                        S_model[*it,*jt] = 0;
                    }
                }
            }
        }
    }
    
    
    
}















void get_cmd(int argc, char* argv[],blockSH & sh,std::string & out)
{
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::varargxxxx<int,int,int,int> iiii(arg,"-iiii","--iiii","iref ia ib ic",1,2,3,4);

        cmd::posarg<std::string> hamla_fn(arg,"input lead hamiltonian file name", "");
        cmd::posarg<int> aoa(arg,"exmol lead atom index overlap", 0);
        cmd::posarg<std::string> hamm_fn(arg,"input exmol hamiltonian file name", "");
        cmd::posarg<int> aob(arg,"exmol lead atom index overlap", 0);
        cmd::posarg<std::string> hamlb_fn(arg,"input lead hamiltonian file name", "");
        
        if(*help)
        {
                arg.print_help(std::cout,"shmat");
                std::cout << "examples:" << std::endl;
                std::cout << "version 0.9b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.9b" << std::endl;
                std::cout << "Copyright (c) 2010 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
        }                                                                                                                


//        if(!~xyz_fn) throw std::runtime_error("Give xyz structure file!");
//        if(!~coup_fn) throw std::runtime_error("Give coupling file!");
//        if(!~iiii) throw std::runtime_error("Give reference atom and lattice neighbours index!");

        std::cout << "ok-1" << std::endl;


        blockSH shla,shlb,shm;
        int ao = *aoa;
        int bo = *aob;
        load(*hamla_fn,shla);

        load(*hamm_fn,shm);
        load(*hamlb_fn,shlb);
    
        for (int i =0; i< shla.orbital_range.size(); i++)
        std::cout << shla.orbital_range[i].first << " " << shla.orbital_range[i].second << std::endl;

        
//        blockSH sh;
        int la = shla.orbital_range.size();
        int lm = shm.orbital_range.size();
        int lb = shlb.orbital_range.size();

        std::cout << "ok0" << std::endl;


        double dea,deb;
        Eigen::MatrixXcd ds,dh,dsm,dhm;
        std::cout << "la:" << la << "    ao:" << ao << std::endl;    
        
        shla.get_block(la-ao,la-ao,ao,ao,ds,dh);
        std::cout << "ok0.1" << std::endl;

        shm.get_block(0,0,ao,ao,dsm,dhm);
        std::cout << "ok0.2" << std::endl;

        dea = std::real((dh-dhm).trace()/ds.trace());

        std::cout << "ok0.5" << std::endl;

        shm.get_block(lm-bo,lm-bo,bo,bo,dsm,dhm);
        shlb.get_block(0,0,bo,bo,ds,dh);
        deb = std::real((dh-dhm).trace()/ds.trace());

        std::cout << "ok1" << std::endl;

        shla.H = shla.H - dea*shla.S;
        shlb.H = shlb.H - deb*shlb.S;
        
        blockSH shmm;
        shm.get_block(ao,la-ao-bo,ao,la-ao-bo,shmm);
        sh = shla;

        std::cout << "ok2" << std::endl;

        sh.add_block(shmm);

        std::cout << "ok3" << std::endl;


        sh.add_block(shlb);

        std::cout << "ok4" << std::endl;

        Eigen::MatrixXcd s,h;

        shm.get_block(0,ao,ao,lm-ao-bo,s,h);
        sh.set_block(la-ao,la,ao,lm-ao-bo,s,h);
        
        shm.get_block(ao,0,lm-ao-bo,ao,s,h);
        sh.set_block(la,la-ao,lm-ao-bo,ao,s,h);

        shm.get_block(ao,lm-bo,lm-ao-bo,ao,s,h);
        sh.set_block(la,la+lm-ao-bo,lm-ao-bo,ao,s,h);

        shm.get_block(lm-bo,ao,ao,lm-ao-bo,s,h);
        sh.set_block(la+lm-ao-bo,la,ao,lm-ao-bo,s,h);

        sh.ef = shm.ef;
}


using namespace std;


int main(int argc, char* argv[])
{
try
{
    blockSH sh;
    std::string outfn;
    get_cmd(argc,argv,sh,outfn);
    
    save(outfn,sh);
    
} catch(const char* message)
{
    std::cerr << message << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
catch(const std::exception & e)
{
    std::cerr << e.what() << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
    return EXIT_SUCCESS;



}




*/

